/* cluster_main.c */

#include "cluster.h"

int main(int argc, char** argv)
{
	int result = 0, return_value = 0, elements_no = 0;
	int g1_size = 0, g2_size = 0;
	int *g = NULL, *g1 = NULL, *g2 = NULL;
	sparse_matrix_lst *A = NULL, *A_g = NULL;
	sub_group *P = NULL, *g_sub_group = NULL, *last_sub_group = NULL;
	sub_group *g1_sub_group = NULL, *g2_sub_group = NULL;
	elem *s = NULL, *k = NULL, *k_g = NULL, *res_vec = NULL;
	elem delta_Q = 0, total_delta_Q = 0;
	elem **B = NULL, **B_g = NULL, **gen_B_g = NULL;
	elem *f = NULL, *empty_f = NULL;
	elem M = 0, percision = 0;
	eigenpair EP;

	static int gcount =0;

	if (argc != 3)
	{
		result = printf("\nusage: cluster <adjacency-mat-file> <precision>\n");
		assert(result >= 0);
		return_value = -1;
		goto main_cleanup;
	}

	A = read_adjacency_matrix(argv[1]);
	if (A == NULL)
	{
		result = printf("\n Could not read the adjacency matrix file \n");
		assert(result >= 0);
		return_value = -1;
		goto main_cleanup;
	}

	percision = atof(argv[2]);

	last_sub_group = create_initial_P(A, &P);
	if (last_sub_group == NULL)
	{
		return_value = -1;
		goto main_cleanup;
	}

	srand(time(NULL));

	B = create_modularity_mat(A);
	k = calc_nodes_degrees(A);
	M = count_edges(k, A->n);

	/*	The moment the first element is indivisible, this means that there are
	no divisible elements  */
	while(P != NULL && P->is_divisible)
	{
		if (last_sub_group == P) {
			last_sub_group = NULL;
		}

		g_sub_group = P;
		P = P->next;

		g = g_sub_group->elements;
		elements_no = g_sub_group->elements_no;
		/* g sub group will reenter p, only if it's indivisible */
		g_sub_group->is_divisible = 0;

		/* Algorithm 2 */
		if(elements_no != A->n) 
		{
			f = create_f(B, g, A->n, elements_no);
			gen_B_g = create_generalized_modularity_matrix(B, g, elements_no,
				f, A->n);
		}
		else 
		{
			/* Optimization for first iteration. */
			f = (elem*)calloc(elements_no, sizeof(elem));
			gen_B_g = B;
		}

		k_g = create_k_g(k, g, elements_no);
		A_g = create_A_g(A, g, elements_no);

		EP = find_eigenpair(A_g, gen_B_g, k_g, f, M, g, elements_no,
			percision);

		if (IS_POSITIVE(EP.eigenvalue))
		{
			int r1 = 0;
			s = (elem*) calloc(elements_no, sizeof(elem));
			assert(s != NULL);

			compute_s(s, g, EP.eigenvector, elements_no);

			res_vec = (elem*)calloc(elements_no, sizeof(elem));
			assert(res_vec != NULL);

			mult_vector_matrix(res_vec, s, gen_B_g, elements_no);
			
			delta_Q = mult_vector_vector(res_vec, s, elements_no) / 2;
			free(res_vec);

			if(elements_no != A->n)
			{
				empty_f = (elem*)calloc(elements_no, sizeof(elem));
				/* Sending zeroed f for B_g and not generalized B_g */
				B_g = create_generalized_modularity_matrix(B, g, elements_no,
					empty_f, A->n);
				free(empty_f);
			}
			else
			{
				/* Optimization for first iteration. */
				B_g = B;
			}

			delta_Q += improve_2_div(B_g, s, k_g, M, elements_no);

			if(elements_no != A->n) 
			{
				free_matrix(B_g, elements_no);
			}

			if (IS_POSITIVE(delta_Q))
			{

				total_delta_Q += delta_Q;
				div_indices(s, g, &g1, &g2, &g1_size, &g2_size, elements_no);
				free(s);

				if (g1_size == elements_no || g1_size == 0)
				{
					add_to_P(&P, &last_sub_group, g_sub_group);
					free(g1);
					free(g2);
				}
				else
				{
					/*	There is no need for the g group, it was divided into
					g1 and g2  */
					free(g_sub_group->elements);
					free(g_sub_group);

					g1_sub_group = (sub_group*)malloc(sizeof(sub_group));
					assert(g1_sub_group != NULL);
					fill_sub_group(g1_sub_group, g1, g1_size);
					add_to_P(&P, &last_sub_group, g1_sub_group);

					g2_sub_group = (sub_group*)malloc(sizeof(sub_group));
					assert(g1_sub_group != NULL);
					fill_sub_group(g2_sub_group, g2, g2_size);
					add_to_P(&P, &last_sub_group, g2_sub_group);
				}
			}
			else
			{
				add_to_P(&P, &last_sub_group, g_sub_group);
			}
		}
		else
		{
			add_to_P(&P, &last_sub_group, g_sub_group);
		}

		if(elements_no == A->n) 
		{
			gen_B_g = NULL;
		}

		if(NULL != gen_B_g) 
		{
			free_matrix(gen_B_g, elements_no);
		}

		free(f);		
		free(EP.eigenvector);
		free_sparse_matrix_lst(A_g);
		free(k_g);
	}

	printf("%lf %d\n", total_delta_Q, count_elements(P));
	print_P(P, A->n);

main_cleanup:
	if (k != NULL)
		free(k);
	if (B != NULL)
		free_matrix(B, A->n);
	if (A != NULL)
		free_sparse_matrix_lst(A);
	if (P != NULL)
		free_P(P);

	return return_value;
}

